Tarea 4: Bayesian Inference

CC6104: Statistical Thinking

Integrantes :

Cuerpo Docente:

Índice:

  1. Objetivo
  2. Instrucciones
  3. Referencias
  4. Elaboración de Código

Objetivo

Bienvenides a la cuarta tarea del curso Statistical Thinking. Esta tarea tiene como objetivo evaluar los contenidos teóricos de la primera parte del curso, los cuales se enfocan principalmente en introducirlos en la estadística bayesiana. Si aún no han visto las clases, se recomienda visitar los enlaces de las referencias.

La tarea consta de una parte práctica con el fin de introducirlos a la programación en R enfocada en el análisis estadístico de datos.

Instrucciones:

Referencias:

Slides de las clases:

Videos de las clases:

Documentación:

#Elaboración de Código

En la siguiente sección deberá resolver cada uno de los experimentos computacionales a través de la programación en R. Para esto se le aconseja que cree funciones en R, ya que le facilitará la ejecución de gran parte de lo solicitado.

Para el desarrollo preste mucha atención en los enunciados, ya que se le solicitará la implementación de métodos sin uso de funciones predefinidas. Por otro lado, Las librerías permitidas para desarrollar de la tarea 4 son las siguientes:

# Manipulación de estructuras
library(tidyverse)
library(dplyr)
library(tidyr)
library(purrr)

# Para realizar plots
library(scatterplot3d)
library(ggplot2)
library(plotly)

# Manipulación de varios plots en una imagen.
library(gridExtra)

# Análisis bayesiano
library("rethinking")

Si no tiene instalada la librería “rethinking”, ejecute las siguientes líneas de código para instalar la librería:

install.packages("rethinking")

En caso de tener problemas al momento de instalar la librería con el código anterior, utilice las siguiente chunk:

install.packages(c("mvtnorm","loo","coda"), repos="https://cloud.r-project.org/",dependencies=TRUE)
options(repos=c(getOption('repos'), rethinking='http://xcelab.net/R'))
install.packages('rethinking',type='source')

Pregunta 1:

Un conjunto de carteros aburridos de las mordidas de perros ha decidido realizar un catastro de mordidas recibidas por los empleados de su empresa en un periodo de 8 meses, planeando en base a estos datos realizar inferencia bayesiana. Los datos de las mordidas estas datos por el dataset no+mordidas.csv, en donde cada fila representa las mordidas recibidas por diferentes carteros y las columnas señalan si fue mordido o no el cartero en los meses de estudio (si fue mordido es señalado con un 1, de lo contrario es señalado con un 0). Cabe señalar que un cartero no puede ser mordido mas de una vez al mes, ya que el damnificado recibe licencia de todo un mes tras la mordida, reincorporándose el siguiente mes al trabajo.

dataMordidas <- read.csv("no+mordidas.csv", header = TRUE)
cols <- c("bites_month_1", "bites_month_2", "bites_month_3",
          "bites_month_4", "bites_month_5", "bites_month_6",
          "bites_month_7", "bites_month_8")

grid_len <- 100
p_grid <- 
prior <- 

posteriors <- c()
counted.months <- c()
p_grids <- rep(p_grid, 4)

for (n in c(1, 2, 4, 8)) {
  x <- # total de mordidas
  
  likelihood <- 
  unstd.posterior <- 
  posterior <- 
  
  posteriors <- c(posteriors, posterior)
  counted.months <- c(counted.months, rep(n, grid_len))
}

post.df <- data.frame("grid"=p_grids, "posterior"=posteriors, "months"=counted.months)

ggplot(post.df) +
  geom_line(aes(x=grid, y=posterior)) +
  facet_grid(months ~ .)

Respuesta:

prior <- 
posteriors <- c()

for (n in c(1, 2, 4, 8)) {
  x <- 
  
  likelihood <- 
  unstd.posterior <- 
  posterior <- 
  
  posteriors <- c(posteriors, posterior)
}

post.df$posterior2 <- posteriors

ggplot(post.df) +
  geom_line(aes(x=grid, y=posterior)) +
  geom_line(aes(x=grid, y=posterior2), color="red") +
  facet_grid(months ~ .)

Respuesta:

for (n in c(1, 2, 4, 8)) {
  res <- quap(
    alist(
      W ~ , # binomial likelihood
      p ~  # uniform pior
    ),
    data = list(W=sum(dataMordidas[, cols[1:n]]),L=n*250-sum(dataMordidas[, cols[1:n]])) )
  
  df <- precis(res)
  curve(dnorm(x, df$mean, df$sd), lty=2, add=FALSE)
}

Respuesta:

¿Cómo puede interpretar los resultados?

grid_len <- 1000
p_grid <- 

prior <- 

likelihood <- 

unstd.posterior <- 

posterior <- 

samples <- sample(p_grid, prob=posterior, size=1e4, replace=TRUE)

dens(samples)

Intervalos:

# fracción de valores <0.35
# fracción de valores entre 0.35 y 0.45
# fracción de valores >0.45

Respuesta:

-[] Calcule un intervalo de credibilidad al \(50%\), \(75%\) y \(95%\). ¿Como se puede interpretar los resultados? ¿Cual podría ser un problema al usar intervalos de credibilidad?

Intervalo \(50\%\):

Intervalo \(75\%\)

Intervalo \(95\%\)

Respuesta:

Intervalo \(50\%\):

Intervalo \(75\%\)

Intervalo \(95\%\)

Respuesta:

ratio.real <- 
samples <- 

ggplot() +
  geom_density(aes(samples)) +
  geom_vline(aes(xintercept=ratio.real), color="red")

Respuesta:

Respuesta:

Pregunta 2:

En esta pregunta se trabajara con el dataset “notas.csv” el cual contiene las notas históricas de un curso desconocido. Suponga que los datos vienen de una distribución \(\mathcal{N}(\mu,\sigma^2)\), el objetivo de la pregunta es estudiar el comportamiento de los datos y los posibles valores de \(\mu,\sigma\) mediante técnicas de inferencia Bayesiana.

Usted sabe un dato extra sobre la información, \(\sigma \in [0.5,1.5]\) y, además, tiene una fuerte creencia a que es mas probable encontrar la desviación estándar real entre \([0.5,1]\) que en \((1,1.5]\).

# Librerias

# Leer información
data_notas <- read.csv("notas.csv")

# Función para crear likelihood dado mu y sigma
grid_function <- function(mu,sigma){
   .... # Funcion de likelihood
}

# Valores de la grilla
grid_mu <- ...
grid_sigma <- ...

# Se crea la grilla 2d
data_grid <- expand_grid(grid_mu,grid_sigma)

# Se guarda la likelihod
data_grid$likelihood <- map2(data_grid$grid_mu,data_grid$grid_sigma, grid_function)

# Se transforma el forma de map2 a una columna
data_grid <- unnest(data_grid,cols = c("likelihood"))


# Valores de los priors
prior_mu <- ...
prior_sigma <-  ...

# Se crea la grilla 2d de priors
prior <- expand_grid(prior_mu,prior_sigma)

# Se calculan los valores del prior
data_grid$prior <-  map2(prior$prior_mu,prior$prior_sigma, prod)
data_grid <- unnest(data_grid,cols = c("prior"))

# Se calcula el posterior
data_grid$unstd_posterior <-  data_grid$likelihood * data_grid$prior

# Se estandariza el posterior
data_grid$posterior <- data_grid$unstd_posterior/sum(data_grid$unstd_posterior)

# Se ajustan los valores de la posterior para que no sean valores tan pequñeos
data_grid$posterior <- (data_grid$posterior - min(data_grid$posterior))/(max(data_grid$posterior)-min(data_grid$posterior))

Respuesta:

# Punto de referencia
# Se recomienda cambiar estos valores por unos adecuados que le permitan estudiar
# Los valores de la distribución de mejor manera
valor_x <- 1
valor_y <- 1

# Grafico

punto_comparacion <- tibble(x = valor_x, y = valor_y)

plt <- data_grid %>%
  ggplot(aes(x = grid_mu, y = grid_sigma)) +
  geom_raster(aes(fill = posterior),
    interpolate = T
  )+ 
  geom_point(x = valor_x, y = valor_y, size = 1.3,color="white")+
  geom_label(
    data = punto_comparacion, aes(x, y),
    label = "Punto Comparación",
    fill = "green",
    color = "black",
    nudge_y = 0, # Este parametro desplaza la caja por el eje y
    nudge_x = 1 # Este parametro desplaza la caja por el eje x
  )+
  scale_fill_viridis_c() +
  labs(
    title = "Posterior para Mean y Standard Deviation",
    x = expression(mu ["Mean"]),
    y = expression(sigma ["Standar Deviation"])
  ) +
  theme(panel.grid = element_blank())

plt

Respuesta:

# Codificamos los datos
x <- 1:length(data_grid$posterior)

# Sampleamos los indices
posterior_samples_aux <- sample(x,size = 1e4, replace = T, prob = data_grid$posterior)

# Obtenemos los verdaderos valores de la sampling distribution
posterior_samples <- data_grid[posterior_samples_aux,]

# Obtenemos solos los valores relevantes para la densidad
df <- data.frame(posterior_samples$grid_mu,posterior_samples$grid_sigma)

# Realizamos las densidades
dens(df)

Respuesta:

Pregunta 3: Regresión Lineal Bayesiana

El objetivo de esta pregunta es introducirlo en la aplicación de una regresión bayesiana. Con esto, buscaremos entender como calcular una regresión bayesiana en base al “motor” aproximación de Laplace, revisando como se comporta la credibilidad de sus predicciones y como la regresión lineal puede llegar a mutar a aplicar una transformación en el vector \(x\). Para responder esta pregunta centre su desarrollo solo en las clases y las funciones que nos brinda la librería rethinking.

Para esta pregunta usaremos el dataset Pokemon.csv, que contiene las características de más de 700 personajes.

poke.df <- read.csv("Pokemon.csv")

Respuesta:

ggplot() +
  geom_point(aes(x=poke.df$Attack, y=poke.df$Sp..Atk), color="red") +
  geom_line(aes(x=, y=), color="orange") +
  geom_ribbon(aes(x=, y=, ymin=, ymax=), alpha=0.3, fill="orange") +
  geom_ribbon(aes(x=, y=, ymin=, ymax=), alpha=0.1, fill="orange")
 

A work by CC6104